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The stochastic dynamics of the damped harmonic oscillator in a heat bath is simulated with 
an algorithm that is exact for time steps of arbitrary size. Exact analytical results are given for 
correlation functions and power spectra in the form they acquire when computed from experimental 
time-lapse recordings. Three applications are discussed: (i) Effects of finite sampling-rate and 
-time, described exactly here, are similar for other stochastic dynamical systems — e.g. motile 
micro-organisms and their time-lapse recorded trajectories, (ii) The same statistics is satisfied by 
any experimental system to the extent it is interpreted as a damped harmonic oscillator at finite 
temperature — such as an AFM cantilever, (iii) Three other models of fundamental interest are 
limiting cases of the damped harmonic oscillator at finite temperature; it consequently bridges 
their differences and describes effects of finite sampling rate and sampling time for these models as 
well. Finally, we give a brief discussion of nondimensionalization. 
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I. INTRODUCTION 

The damped harmonic oscillator in a heat bath is 
the archetypical bounded Brownian dynamical system 
with inertia and the simplest possible of this kind. An- 
alytically solvable, it offers insights that are valid also 
for more complex systems. Its experimental correlation 
functions and power spectra are given analytically here 
in the form they take when computed from time-lapse 
recorded trajectories. Such trajectories are generated 
and analyzed for illustration in an exact Monte Carlo 
simulation. The results differ significantly from those 
derived from the standard continuous-sampling formu- 
lation. 

In mathematical terms, it is the Ornstein-Uhlenbeck 
(OU) process in a Hookean force field we treat, and the 
results discussed here extend the ones given in [19, 23] as 
well as the free case treated by Gillespie [11]. Three clas- 
sical papers on the OU-process are [7, 29, 30]. Histori- 
cally, the model was prompted by a question by Smolu- 
chowski regarding how inertia might modify Einstein's 
theory of Brownian motion [20]. H. A. Lorcntz soon af- 
ter pointed out that Ornstein's answer to that question 
was insufficient for the case of classical Brownian mo- 
tion, i.e., in a fluid of density similar to the Brownian 
particle [16]: Hydrodynamical effects, entrainment and 
back-flow, were more important with the time resolu- 
tion that was available in the first century of Brownian 
motion. Inertia in classical Brownian motion became 
experimentally relevant only with precision calibration 
of optical tweezers [2, 3] and was directly observed only 
in 2005 [17, 27]. 

The OU process remained and remains, however, an 
important model for Brownian motion dominated by in- 
ertia, such as massive particles [5, 6, 15] or AFM can- 
tilevers [22] in air, and for other kinds of persistent ran- 



dom motion, e.g. cell migration [8-10, 12-14, 25, 26] and 
commodity pricing [24]. Consequently, the treatment 
given here of its experimental statistics for time-lapse 
recorded data should be useful in several ways: 

(i) All experimental statistics contain effects of finite 
sampling rate, finite sampling time, and finite statis- 
tics. They do so to various degrees, but the effects are 
inherent in measurements. We give exact analytical ex- 
pressions for these effects, for the model treated here. 
Statistics for more complex systems contain the same 
effects, qualitatively, and quantitatively as well, by de- 
grees we estimate. 

(ii) The statistics we describe below must be found 
for any experimental system to the extent that system 
is interpreted as a damped harmonic oscillator at finite 
temperature — e.g. an AFM cantilever to be calibrated 
by interpretation of its thermal power spectrum. 

(iii) Three other models are limiting cases of the 
model treated here: 

1. At vanishing mass, Einstein's theory for Brownian 
motion in a harmonic trap, which, e.g., is a min- 
imal model for the Brownian dynamics of a mi- 
crosphere held in an optical trap, with magnetic 
tweezers [28], or surface-tethered by DNA [1]. 

2. At vanishing external force, the Ornstein- 
Uhlenbeck model of free Brownian motion with 
inertia, which, e.g., is a minimal model for the 
persistent random motion seen in trajectories of 
motile cells. 

3. At vanishing mass and external force, Einstein's 
original theory for Brownian motion in a fluid at 
rest. 

The results for the harmonic oscillator, given below, 
carry over to these three models. The limits are not 
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all obvious, but always enlightening, hence described 
below. 

(iv) As the model treated here bridges the differences 
between the three limiting cases, the material pre- 
sented here is well suited for a pedagogical, hands-on 
computer-based introduction to the four dynamic sys- 
tems covered here: free/bound diffusion, with/without 
inertia, their equilibrium behavior, correlations, and 
power spectra, and their transient behavior to equilib- 
rium. 



II. EXACT DISCRETIZED 
EINSTEIN-ORNSTEIN-UHLENBECK THEORY 
OF BROWNIAN MOTION IN HOOKEAN 
FORCE FIELD 

The Einstein-Ornstein-Uhlenbeck theory for the 
Brownian motion of a damped harmonic oscillator in one 
dimension is simply Newton's Second Law for the oscil- 
lator with a thermal driving force, a.k.a. the Langevin 
equation for this system, 

mx(t) + J±(t) + Kx(t) = Fthcrm(i) (1) 

Fthcrm(t) = (2fc B T 7 ) 1 / 2 ry(t) . (2) 

Here x(t) is the coordinate of the oscillator as function 
of time t, m its inertial mass, 7 its friction coefficient, n 
is Hooke's constant, and F t horm is the thermal force on 
the oscillator. Equation (2) gives the amplitude of this 
thermal noise explicitly in terms of 7, the Boltzmann 
energy k B T, and r](t), which is a normalized white-noise 
process, i.e., the time derivative of a Wiener process, 
77 = dW/dt, hence 

(j7(i))=0; (v(t)v(t)) =6{t-t') for all t,* 7 . (3) 

Equation (1) can be rewritten as two coupled first-order 
differential equations, 

dt\v(t))- M {v(t))+ T {ri(t)) ' ^ 

where we have introduced Einstein's relation D = 
knT/^i and the 2x2 matrix 

Here ujq = yjn/m is the cyclic frequency of the un- 
damped oscillator, and r = m/ 7 is the characteris- 
tic time of the exponential decrease with time that 
the momentum of the particle undergoes in the ab- 
sence of all but friction forces. Below, we shall also 
need the cyclic frequency of the damp ed oscillator, 
uj = njm — 7 2 /(4m 2 ) = \/uJq — 1/(4t 2 ), which is real 
for less than critical damping, 7 2 < 4m«. 
Equation (4) is solved by 

(:|iO^/> M, '-° U)- < 6 » 



which, for arbitrary positive At, and with tj = jAt, 
xj = x(tj), and Vj — v(tj), gives us the recursive relation 

(?;;)=- Ma *(S ■ m 

where 

(8) 

and the time-independent matrix exponential can be 
written 

e -MAt =e -% [cos^A^j + sin ( w At) j] (9) 

with 

is G?) andJs (S i) • (10) 

That the matrix exponential can be written this way, 
can be proven in several ways. We used the algebra 
of Pauli matrices. Alternatively, one may observe that 
the cosine and sine terms on the RHS are the even and 
odd parts of the LHS. The latter are straight-forwardly, 
if tediously, computed from the Taylor series for the 
exponential. Inserting Eq. (9) in Eq. (8) we see that 

V2D f ti+1 tj+i-t 
Axj = J Ate 2- sm(u(tj +1 - t)) r)(t) (11) 



V2D r i+1 tj+i-t 
Av j = -7^2 J dte 2T sm(u(t j+1 - t)) f](t) 

(12) 

V2D f ti+1 tj+i-t 
H dte 2, cos(w(tj + i — t)) t](t) 

are two correlated random numbers from zero-mean 
Gaussian distributions. They can be written as a lin- 
ear combination of two independent Gaussian variables: 
The four parameters that determine this linear combi- 
nation can be chosen at will, as long as the combination 
has the same variance-covariance as the two original cor- 
related variables. This is a direct consequence of the 
Gaussian distribution being completely determined by 
its mean and variance-covariance. Thus, we can write 

Axj = a xx i", (13) 
A «j = <>lv/°xx (.j + \Ztfv-<rtv/°L tj . ( 14 ) 

where the cts are elements of the variance-covariance ma- 
trix (see below), and £ and £ are independent random 
numbers with Gaussian distribution, unit variance, and 
zero mean. This particular choice of linear combina- 
tion mirrors the structure of Eqs. (11) and (12). Using 
Eqs. (3), (11), and (12) we calculate that the elements 



of the variance-covariance matrix are, for uj 7^ 
D 



a 2 xx = ((Axj) 2 ) 



(4w 



2^2 



(15) 



e ^ [cos(2wAi) - 2ujt sin(2wAi) - 4w 2 r 2 " 
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al ee ((Avj) 2 ) = ^ (4 W 2 r 2 + 



(16) 



e ^ [cos(2wAi) + 2wr sin(2wAi) - 4lu 2 t 2 ' 



°lv = (^XjAvj) = -=-je ? sin 2 (wAi) 



(17) 



At critical damping (ui 
simplifies to 

al x = ((Axj) 2 ) 



0) this variance-covariance 
(18) 



1 



= ADt 1 - e~ [1 + At It + -(At/r) 2 } 



((Av 3 f) 

D/t (l - [1 - At/r + -(At/r) 2 ] 



,2 _ 



AxjAvj) = De 



(At/r) 2 



(19) 



(20) 



Figure 1 shows the simulated positions in the un- 
derdamped, critically damped, and overdamped regime. 
Note that all three regimes have the same mean (zero) 
and variance (k^T/n). Only at short times, as shown 
in the insert, does the difference between the three regi- 
mens reveal itself: The position coordinate of the under- 
damped system oscillates, while the position coordinate 
of the overdamped system does not show discernible per- 
sistence of motion. The position coordinate of the criti- 
cally damped system looks at times like it oscillates, but 
actually only displays positively correlated random mo- 
tion, as revealed by the velocity auto-correlation func- 
tion given below and shown in Fig. 3B. 



III. MEAN-SQUARED-DISPLACEMENT 

When a time-series of positions is examined, the 
first measure applied is frequently the mean-squared- 
displacent (MSD). This has to do with its ease of com- 
putation, the existence of exact analytical results for 
the expectation value, and comparative robustness of 
the measure to experimental measurement errors. Also, 
there are no discretization effects to worry about — the 
MSD for finite sampling frequency is simply the MSD for 
continuous recording, evaluated at discrete time points. 
More importantly, unlike the covariance and the power- 
spectra discussed below, the MSD is well-defined even 
if the process is not bounded. In the present example, 
both the position and velocity processes are bounded, 
and we therefore have a simple relation between the 
MSD and the auto-covariance: 

MSD(i) ee ((x{t) - x(0)) 2 ) = 2(x 2 ) - 2<a;(t)x(0)> (21) 
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FIG. 1. Positions of harmonically trapped massive parti- 
cles in a thermal bath. Three different regimes are studied 
by varying the drag coefficient above and below the crit- 
ical value, where to = 0, by an order of magnitude: Un- 
derdamped (red, 7 = 7 crlt /10), critically damped (green, 
7 = 7 crlt , offset by 2nm), and overdamped (blue, 7 = 
107 crl , offset by 4nm). Simulation parameters in Eq. (7): 
m = lng, K = 225mg/s 2 , T = 275 K, / samp i c = 65, 536 Hz, 
(x(0),v(0)) = (0,0), and the same (£j,r)j) are used for the 
three simulations. The values for the underdamped case are 
alike to those for an AFM cantilever in water, whereas the 
critical and overdamped case corresponds to the same can- 
tilever in a more viscous environment or with a smaller spring 
constant. The insert shows a magnified portion of the trace, 
revealing the oscillating, critical, and random nature of the 
motion, respectively. On the right, histograms show the dis- 
tribution of the position data with &b T/k- variance Gaus- 
sians overlaid. 



re [x") — k^T/n, i.e., the MSD is a constant minus 
twice the auto-covariance of the position process (see 
Sec. tV). This is illustrated in Fig. 2 where we show the 
MSD for the three time series plotted in Fig. 1. 



IV. COVARIANCE 



The stationary-state solution given in Eq. (6) oscil- 
lates in resonating response to the thermal noise that 
drives it, which is seen by its oscillating auto-covariance 
functions, but its oscillations are randomly phase shifted 
by the same noise that drives oscillations, which causes 
the exponential decay in auto-covariance functions, see 
below. The result is that equal-time ensemble averages 
like (x(t) 2 } are constant in time. This is easily proven 
for quadratic expressions, either directly from the solu- 
tion given in Eq. (6), or by differentiation w.r.t. time 
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FIG. 2. Mean-squared-displacements for harmonically 
trapped massive particles. Same simulation conditions as 
in Fig. 1, except t msr — 8 s. The result of the stochastic 
simulations are shown in red, green, and blue. Comple- 
mentary colors (cyan, magenta, and yellow) show the MSD, 
Eq. (21), evaluated at discrete time-points with the same pa- 
rameter values as in the simulation. Thin black lines show 
the MSD Eq. (21) for continuous time. Red and blue lines 
show slopes of two and one, respectively. Notice how the 
MSDs all plateau at the same level for time-lags larger than 
3 s. Only for shorter time-lags do the three regimes differ. 



using Ito's lemma: 
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As is seen by inspection, this equation has the time- 
independent solution 



(23) 



in accord with the equipartition theorem. This solution 
is the unique attractor for the system's dynamics: Left 
to itself, any discrepancy from this time-independent so- 
lution will decrease to zero exponentially fast in time — 
possibly while oscillating harmonically with cyclic fre- 
quency 2u> — as is seen from the fact that the 3 x 3-matrix 
in Eq. (23) has eigenvalues — 1/r and — 1/r ± i2cj, and 
determinant — Auq/t. 

Similar reasoning (or simply multiplying Eq. (7) with 
(xj, Vj), then taking the expectation value, and applying 
the equipartition theorem) gives the covariances 



(x(t)x(0)) (x(t)v(0)) 
(v(t)x(0)) (v(t)v(0)) 



(24) 



D 



cos Lji+ - 



For uj real (underdamped system), these covariances os- 
cillate, we see, with amplitudes that decrease exponen- 
tially in time with characteristic time 2r. For u) imagi- 
nary (overdamped system), we rewrite Eq. (24) as 



(x(t)x(Q)) (x(t)v(Q)} 

M*M0)> (v(t)v(o)) 



(25) 



D 

2Ur 



4 T + 



-t/T_ 



1 =i 



1 -i 



from which we see that the system decreases as a double- 
exponential with characteristic times t± = 2t/(1 ± 
2r|u;|). At critical damping (w = 0), the expressions 
for the covariances simplify to 



(x(t)x(0)) (x(t)v(0)) 
(v(t)x(0)) <«(i)«(0)> 



(26) 



D 



l+t/(2r 
-t 



1 



t 

t/(2r) 



Figure 3 illustrates these dynamics for the normalized 
covariances, aka the correlation functions. 

As is the case for the MSD, there are no discretiza- 
tion effects to worry about — the covariance for finite 
sampling frequency is simply the covariance for con- 
tinuous recording, evaluated at discrete time points. 
However, the given velocity correlations, and thus also 
the position-velocity correlations, are only correct for 
the actual, instantaneous velocities: If the instanta- 
neous velocity is not directly measured, but instead es- 
timated from measured positions as a "secant-velocity" 
(see Eq. (60) and Sec. VII A 1), then the correspond- 
ing secant-velocity correlation function should be cal- 
culated from the position covariance function given in 
Eq. (24). Figure 4 shows how such secant- velocity cor- 
relations deviate from the instantaneous-velocity corre- 
lations, mainly for short time-lags. 



V. POWER SPECTRA 

The power spectral density (PSD) of a variable x(t) 
or v(t) is defined as the expectation value of the squared 
modulus of its Fourier transform, with a normalization 
that differs between authors. We choose one in which 
the power spectral density remains constant for increas- 
ing measurement time, t msl - 



and 



p(")(/ fc ) = (\v k \ 2 )/t msl = (27r/ fe ) 2 P^(/ fe ) 
2D 



(2ttt)2(A--/o7/) 2 + 1 ' 



(28) 



COS LUt — 



2ljt 



where 2tt/q — ujq is the resonance frequency of the 
system. Here, we used Vk = i2nxk, which is an ap- 
proximation that ignores contributions from the finite 
measurement-time. 
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FIG. 3. Correlation functions, Eq. (24), for massive par- 
ticles in a harmonic potential, driven by thermal noise. 
Color-coding and simulation settings are the same as in 
Fig. 2: Red, green, and blue show the result of a stochas- 
tic simulation; complementary colors show the covariance 
(not a fit), Eq. (24), normalized and evaluated at dis- 
crete time-points; thin black lines show the covariance for 
continuous time (not a fit); and dashed red lines show 
the enveloping exponential exp(— t/(2r))). Panel A: Po- 
sition auto-correlations, {x(t)x(0)} / {x 2 ). Panel B: Veloc- 
ity auto-correlations, (v(t)v(O)} / {v 2 } . Panel C: Position- 
velocity cross-correlations, (x(t)v(0)) / \J (x 2 )(v 2 ). The 
velocity-position cross-correlation, (v(t)x(0)) / \J (x 2 )(v 2 ), 
(not shown) has the opposite sign but is otherwise identi- 
cal. 
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FIG. 4. Secant- velocity autocorrelation functions for massive 
particles in a harmonic potential, driven by thermal noise. 
Thin black lines, red dashed lines, and simulation settings 
are the same as in Fig. 3B. Additional symbols: Cyan-red, 
magenta-green, and yellow-blue squares show the autocorre- 
lations for Wj — (xj — Xj~i)/At calculated using (x(t)x(0)) 
from Eq. (24). Insert shows the first 0.1 ms on a linear time- 
scale, revealing the very fast decay of the overdamped os- 
cillator and how the secant-velocity autocorrelations over- 
estimate the instantaneous-velocity autocorrelation decay- 
times. As expected, the secant-velocity does a better job 
estimating the instantaneous velocity for the less damped 
system because of its smoother particle-trajectory. 



The finite-time Fourier transform (FTFT) used above 
is denned as 



dte i2nfkt x(t), f k = k/t msl , k integer 



(29) 

In a similar manner we can calculate the finite-sampling- 
frequency PSDs by discrete Fourier transformation 
(DFT) of Eq. (7) 



T/,//N yq I x k \ = g— MA* ( Xk ) + f 

Vk ' \ Vk \ Aw fe 



(30) 



where 



N 



z k = AtJ2e l2 ^ k/N Zj = ^_£V*iA/A^ , (31) 



N 



3 = 1 

for z = x, v, Ax, Av; xn+i — x\ and un+i = i>i; 

fk = kAf = k/t msl 
fc = 0,l,2,...,iV-l 

^msr N At — N //sample 
^ * - /sample 



/Nyq = ~^ A f = 



(32) 

(33) 
(34) 

(35) 



and, 



(Ax* k Ax k ,) = a 2 xx (t k ik>) = °l x At W h# (36) 



(Av* k Av k ,) = (jy V Att msi S k ,k' 
(Ax* k Av k ,) = a 2 xv Att msr 5 k: k> 



(37) 
(38) 



G 



After isolating (x k ,v k ) in Eq. (30), we find 



(\i k \ 2 ) 



2 

® ' XV 



(A 



where the 2x2 matrix and its inverse are 



a k I + (33 : 



A- 1 



Ofel - (33 



with 



a k = e ~ infk/f ^ - cos(wAt) 
(3 = — sin(wAi) , 



(39) 
(40) 



(41) 
(42) 



complex scalars ((3/uj is real), Hermitian conjugation is 
denoted by f , and complex conjugation by *. 

For the discrete positional power spectrum we thus 
get (see also [19, Appendix A]) 



P ( k x) = <M 2 )A 

\a k 



_1_\2 2 
2ujt i xx 



2^-Re\a k 



(43) 



' XV 1 uj 2 vv 



K +/?T /sample 

whereas the discrete velocity power spectrum is 

>(") — /!„-._ |2\ 



p?> = (\v k n/tn 



(44) 



\a k + ^| 2 ^ + 2^Rc{ Qfc + ^}a 2 xv + 



,4 f) 



sample 

In case of critical damping (w = 0) we simply in- 
sert /3/ui = —e~^At in Eqs. (43) and (44), and use 
Eqs. (18,19,20) for the variance-covariances. 

Figure 5 shows the power spectra that result from nu- 
merical simulations of AFM cantilever positions and ve- 
locities using Eq. (7) for three different drag coefficients, 
as well as the corresponding analytical expressions as 
given in Eqs. (43,44). 

An alternative route to these discrete PSDs is to take 
the discrete Fourier transform of the covariance function 
from the previous section: Both the position and the 
velocity processes are stationary (the joint probability 
density function for each is independent of time) , so the 
Wicncr-Khinchin theorem applies and the Fourier trans- 
form of the correlation function is equal to the PSD. For 
vanishing spring constant, k, the position process is un- 
bounded (not stationary) and the covariance, as well 
as the PSD, are consequently ill-defined. This limit is 
treated in Sec. VII A. 

Note, the PSD for discretely measured (instanta- 
neous) velocities, Eq. (44), should not be confused with 
the discrete PSD for secant- velocities, Eq. (61). The 
latter is the correct expression to use in experiments 
where the velocities are estimated from the measured 
positions. 



VI. VANISHING At: THE LIMIT OF 
CONTINUOUS RECORDING 



As At — > we find to first order in At that the co- 
variances in Eqs. (15-20) reduce to = 0, cr^v — 
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FIG. 5. Power spectra for positions and velocities of harmon- 
ically trapped massive particles in thermal bath, normalized 
by 7. The color-scheme is the same as in Figs. 2 and 3: 
Synthetic data is shown in red, green, and blue; the aliased 
theory, Eqs. (43) and (44), is shown in complementary col- 
ors; and the non-aliased theory, Eqs. (27) and (28), as thin 
black lines. Simulation parameters are as in Fig. 1 (giving 
luo — 15 kHz) with n w in = 32 Hann windows applied to the 
imsr = 8s long time-series. Panel A: Power spectra for po- 
sitions, f~ 2 and /~ 4 behaviors are illustrated by the thin 
blue and red lines, respectively. Notice the disagreement 
between the non-aliased (continuous recording) theory and 
data at high frequencies. Panel B: Power spectra for veloc- 
ities. Dashed black lines show aliased versions of Eq. (28), 
i.e., P<«.")(/) = E^-oo-P W (/ + n/ samplc ), [2], here 
truncated to 201 terms. The non-aliased theory severely 
underestimates the power at high and low frequencies, and 
completely misses the data in the overdamped case. 



and 



= 0. That is, Axj = and 



2DAt/i 

Avj — UwCj — s/lDAtjr Q. In other words, the ve- 
locity process is seen to be driving the position process. 

In this limit the positional PSD takes on the famil- 
iar form given in Eq. (27), and the velocity PSD is 
given in Eq. (28). Comparing the continuous record- 
ing PSDs Eqs. (27,28) with their discrete sampling ana- 
logues Eqs. (43,44) we see in Fig. 5 A that they differ 
substantially at high frequencies, and in Fig. 5B that 
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they differ everywhere but near / in the underdamped 
case. Thus, one should not attempt to fit a theoretical 
PSD derived for "continuously recorded" trajectories to 
an experimental PSD, which necessarily is obtained with 
time-lapse recording. More specifically, one should only 
attempt to fit it to the low-frequency part of the experi- 
mental spectrum for positions, or account for "aliasing" 
before fitting, as described e.g. in [2, Appendix H], which 
turns Eqs. (27) and (28) into Eqs. (43), respectively (44). 

The position-covariance and the MSD do not suffer 
from this dichotomy. None of them change their shape 
or form when switching between discrete and continuous 
time: The discrete versions are equal to the continuous 
versions, evaluated at discrete times: 



(xjXj+t) = (x(tj)x(t j+ e)) = (x(0)x(ti)) 



(45) 



and 



((xj-Xj+t) 2 ) = ((x(tj)-x(t j+t )) 2 ) = ((x(0)-x(ti)) 2 ) , 

(46) 

where, in the last step, we used that the process is sta- 
tionary and the measures therefore invariant to time- 
translations. That is, these measures are unaffected by 
the unavoidable discreteness of real- world data. 

The velocity's covariance is also unaffected by time- 
lapse recording, if we can measure the instantaneous 
velocity, i.e. the velocity vector that is tangential to 
the trajectory of the position. If we cannot and time- 
lapse record only positions, we have a different situation, 
which is treated below. 



VII. VARIOUS PHYSICAL LIMITS: A 
REFERENCE-SET OF FORMULAS 

Three physical limiting cases are of particular inter- 
est: (i) Vanishing spring constant, k = 0, which is the 
Ornstein-Uhlcnbeck theory for Brownian motion for a 
free particle with inertia; (ii) vanishing mass, m = 0, 
which is Einstein's theory for Brownian motion of a par- 
ticle trapped by a Hookean force, and a popular mini- 
malist model for the Brownian motion of a microsphere 
in an optical trap; and (iii) vanishing mass and spring 
constant, which is Einstein's theory for Brownian mo- 
tion of a free particle. 



A. Vanishing spring constant: 
The Ornstein-Uhlenbeck process 

When there is no Hookean restoring force k = 0, 
Eq. (1) describes the free diffusion of a massive particle 
according to Ornstein and Uhlenbeck [29] 



mv(t) + jv(t) = F t herm(i) 



(47) 



This velocity-process is known as the Ornstein- 
Uhlcnbeck (OU) process [29]. When modeling other 
dynamical systems, such as migrating cells, m does not 
refer to the physical mass of the cell but rather its inertia 
to velocity-changes, or persistence of motion; likewise 7 



is not the friction between the cell and substrate but de- 
scribes the rate of memory-loss for the velocity process. 
The structure of Eqs. (7,9,10,13,14) remain the same, 
the only change is uj — in Eq. (10), hence to = i/(2r), 
and the covariances, Eqs. (15,16,17), consequently re- 
duce to 



XX 

2 



Dt (2At/r - 4(1 
D/t (1 - a 2 ) 
D(l-af , 



a) + (l- a 2 )) 



(48) 
(49) 
(50) 



where 



,-At/r 



m/7 . 



(51) 



For the velocity and position processes we thus find 



v j+1 = avj + Avj 



= Xj + t(1 — a)Vj + Aa 



(52) 
(53) 



That is, the velocity process has reduced to a stable 
autoregressive model of order one, AR(l). Its time in- 
tegral, the position process, is unbounded — the particle 
is diffusing freely and without limits. Exact numeri- 
cal update formulas for Vj and Xj have previously been 
given in [11], here we re-derived them for consistency of 
notation. 

(x) 

The discrete-time positional PSD, P k , is obtain- 
able from Eq. (43) with u> — and the as given in 
Eqs. (48,49,50); the expression does not simplify sig- 
nificantly compared to Eq. (43). The discrete-time 

(v) 

velocity PSD, P k , is much simpler and can be de- 
rived directly from the discrete-time velocity process, 
Eq. (52); or by taking the k = limit of Eq. (44). From 
these discrete-time PSDs the continuous-recording ex- 
pressions, P( x \f) and P( v \f), can be obtained by ex- 
panding to leading order in At/r and k/N — f k //sample! 
or they can be derived directly from the continuous-time 
equation of motion Eq. (47) by Fourier transformation 
and, for the positional PSD, remembering that x = v. 
The continuous-recording positional PSD and the two 
velocity PSDs then read 



P {x \h) 
P (v) (fk) = 



D/(2n 2 ) 



(V) _ Uyy Af 



l + a 2 -2acos(27rfc/A0 
2D 



i + (2TTf kT y ' 



(54) 
(55) 
(56) 



Note, that the positional PSD diverges in f k = due 
to the process being unbounded. This implies that we 
cannot use the Wiener-Khinchin theorem to obtain the 
covariance by Fourier transforming the PSDs, and vice 
versa. The velocity process is bounded, however, and 
the velocity correlation functions are straightforward to 
calculate, so we simply list the results for the discrete 
and the continuous case 



(ViVj) 



D 

T 

D 



(v(t)v(t>)) = - e 



-\i-j\At/r 



-\t-t'\/r 



(57) 
(58) 
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As already mentioned, the time integral of the OU- 
process is one of the instances where the mean-squared- 
displacement provides a useful measure for the position 
process, while the positional covariance function is ill- 
defined. The mean-squared-displacement of the time 
integral of the OU-process is 



((x(t) - x(0)) 2 ) = 2Dt U/t + e- t/r - l) 



(59) 



For i<r this MSD increases as t 2 (ballistically) and for 
t 3> t as t (diffusively), with an exponential cross-over 
between the two regimes with characteristic time r. 



1. Secant-velocities 

In experiments, the velocity is typically not measured 
directly, but approximated from the measured positions 
as a "secant- velocity" : 



Cj-l 



)/At 



(60) 



P 



(,:) 



(61) 



with discrete power spectral density 

pW = (\w k \ 2 )/t msi 

= 2(1 - cos(7r/ fc // Nyq )) 
(At) 2 

This is a general result that follows from the definition 
of Wj and is independent of the dynamic model. 

For the OU-process, we can relate this PSD to the 
PSD for the continuous velocity that it approximates, if 
we rewrite it as 

Pl w) =P ( k v) (r/At) 2 (l-a) 2 +a 2 xx /At 

+ r(l - a) Ax k ) + c.c.) . (62) 

Here the first term is proportional to pW, the second 
term is a constant, and the last term is frequency depen- 
dent with c.c. the complex conjugate of the other term 
in the bracket. However, it is not easy to see from this 
expression how much pw deviates from p( v \ We can 
get a good idea about the shape of the PSD from the 
auto-correlation function 



2(v 



1 - \At/r 



(At/r - 1 



and for I > 



(^) (cosh(At/r) 



1 



(63) 
(64) 

1) (65) 
(66) 



where we used Eqs. (13,14,50,53,57) to derive 
Eqs. (63,65) — the latter two expressions are propor- 
tional to the velocity autocorrelations they approximate, 
they only differ by multiplicative factors that are inde- 
pendent of the time-lag for I > 0. To first order in At/r 
we thus have, for £ > 0, 



10 



10 



o 

Q_ 
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OU secant velocity, simulation 
OU velocity, theory 
OU secant velocity, theory 
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FIG. 6. Power spectra of velocities and secant-velocities for 
the OU process. Blue: Simulated velocities, Eq. (52). Red: 
Secant-velocities, Eq. (60), calculated from simulated posi- 
tions, Eq. (53). Yellow: Aliased velocity theory, Eq. (55). 
Cyan: Aliased secant- velocity theory, Eq. (61). Black line: 
Approximated aliased secant-velocity theory, Eq. (68). Sim- 
ulation settings are given the caption of Fig. 7. Notice the 
discrepancy between the true velocity and the estimated (se- 
cant) velocity at high frequencies. 



We now apply the Wiener-Khinchin theorem to get the 
PSD as the Fourier transform of this approximated auto- 
correlation function and find 



,(f) 



P, 



O) 



1 f At 

3 



D 



This approximation to the secant-velocity PSD, as well 
as the exact expression Eq. (61), are both shown in 
Fig. 6, together with numerical simulation results and 
the PSD of the continuous velocity Eq. (55). Notice, 
that even with At/r = 0.27 < 1 (At = l// saro ple and 
r = m/7), the relative difference between P^ and P^ 
is substantial for higher frequencies. Thus, proper care 
should be taken if empirical testing of a model includes 
the fitting of a theoretical velocity PSD to an experimen- 
tal secant-velocity PSD, even when aliasing is properly 
accounted for. But, we now also know how to do this 
correctly: Either add l/3(At/r) 2 D to the experimental 

secant-velocity PSD before fitting with pf? , or fit di- 
rectly using the full theoretical expression for Pl w) given 
in Eq. (61). In cases where this corrective constant is 
much smaller than the velocity PSD for all frequencies 
fitted, it can obviously be ignored. 

We observe that the secant velocity is the time- 
average of the real velocity in the time interval spanned 
by the secant. Time-averaging is low-pass filtering, as is 
well known [21, Chapter 13] and demonstrated in Fig. 6. 



1 lAt X 
3 r 



(67) 



B. Vanishing mass: 
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The optical trap limit 



hence 



The position process in the limit of vanishing mass 
is mathematically identical to the OU velocity process 
in the limit of vanishing trapping force treated above. 
Equation (1) reduces to 



(0) = o, fo^) = hi 



(81) 



jx(t) + Kx(t) = F thcim (t) 



with solution 



!(*) = - f dt' e - 2 -/c(*-t') Ftherm(t ') 

7 J -co 



(69) 



(70) 



where the corner frequency / c = k/(2itj) is the fre- 
quency where P {x) {f) = P^(0)/2, see below. Recy- 
cling the results from Section VII A we can directly write 
down the autocovariance 



(x(t)x(t')) = (x 2 ) e - 2 -/d*-*'l 
with (x 2 ) = k^T / k, and the discrete update rules 

Xj + i = exj + Axj 
c = cxp(— nAt/j) 



(71) 



(72) 
(73) 

(74) 



where £ are uncorrelated random numbers of unit vari- 
ance, zero mean, and Guassian distribution. Likewise, 
the position PSDs for discrete sampling, respectively, 
continuous recording are found to be 



>(*) 



P, 

P [x) Uk) = 



D(l- c 2 ) Atj/n 
' ~ 1 + c 2 -2ccos{2Trk/N) 
D/(2n 2 ) 



fc+f k 



(75) 
(76) 



Finally, the mean-squared-displacemcnt is 

((x(t) - x(0)) 2 ) = 2Dj/k (l - e-*T/«) , (77) 

which shows an exponential cross-over from linear de- 
pendence on t to the constant value 2Dj/k as t >• n/-f, 
see Fig. 2. 



C. Vanishing mass and spring constant: 
Einstein's Brownian motion 



and the velocity and position PSDs for discrete and con- 
tinuous sampling take on the simple forms: 



pO) _ 



P {x) (fk) = 



-^/./sample 



2 sin (7r/ fe // samplo ) 
D 



pjy) =P (v) {fk) = 2D 



(82) 

(83) 
(84) 



As was the case for the time integral of the OU process, 
the position PSDs are singular in fk — because the 
process is unbounded. That is, the meaningful measure 
to study is not the covariance, but rather the mean- 
squarcd-displacement 



{(x(t) - x(0)) 2 ) = 2Dt 



(85) 



which is one of the few well-defined statistics for Ein- 
stein's theory for Brownian motion: Because the trajec- 
tory of positions is a fractal, attempt at estimating the 
average speed of Brownian motion from the displace- 
ment of position occurring in a given time-interval will 
depend on the duration, At, of this interval as l/y/A~t, 
hence diverge when accuracy is sought improved by re- 
ducing At. This was not appreciated before Einstein's 
1905 paper on the subject. 

Figures 7 and 8 show the power spectra, respectively 
mean-squared-displacements obtained in numerical sim- 
ulations of free diffusion and trapped diffusion, as well 
as the graphs of the corresponding analytical expres- 
sions. At short time-scales, i.e., at high frequencies 
in Fig. 7 and for small time-lags in Fig. 8, the ther- 
mal forces dominate and the Hookean force has not had 
time to influence the motion through its constant, but 
weak, confining effect. That is why the Brownian mo- 
tion (green) and optical trap (red) data collapse in this 
regime, whereas the Ornstein-Uhlenbeck process (blue) 
differ from the two due to inertial effects. Conversely, 
at long time scales, low frequencies in Fig. 7 and large 
time-lags in Fig. 8, inertia plays no role, so the Ornstein- 
Uhlenbeck process and Einstein's theory for Brownian 
motion are indistinguishable, whereas the Hookean force 
has had time to exert its confining effect on the optical 
trap data. 



When both m = and k = 0, Eq. (1) reduces to 

yx{t) = F thcrm (t) (78) 

so that 



where 



= xj + V2DA~i £j 



I rtj+i 

VaI L 



(79) 



(80) 



VIII. NONDIMENSIONALIZATION 

Above, we worked with dimension-full equations be- 
cause they are physically intuitive and allow for dimen- 
sional checks of the calculations along the way. Extra 
insight can be gained however, by nondimensionalizing 
the equations. 

As an example, in the dynamic equation for the OT 
case 



7± + kx = F thcim (t) = ^2k B T-fr){t) 



(86) 
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io 1 10 2 10 3 10 4 

Frequency (Hz) 

FIG. 7. Power spectra of positions. Green points: Freely dif- 
fusing massless particle (Einstein's Brownian motion); Red 
points: Trapped massless particle (OT limit, or OU veloc- 
ity process); and Blue points: Freely diffusing massive par- 
ticles (time integral of OU process). Complementary col- 
ors show the finite sampling frequency (aliased) theories, 
whereas the continuous recording (non-aliased) theories are 
shown as thin black lines. Green and blue lines indicate 
f~ 2 and f~ 4 behavior, respectively. Simulation parameters: 
D = 0.46 fim 2 /s, T = 275 K, m = 1 ng, and f c = 500 Hz, 
/sample = 32, 768 Hz, N = 131,072, n win = 32 Hann win- 
dows. The simulation parameters for the OT case are those 
of a 1 fim diameter polystyrene sphere held in an optical trap 
in water at room temperature. The parameters for Einstein's 
Brownian motion are the same, except n = 0. For the OU 
process we increased the density of the sphere roughly 2,000 
times, which is not a physically realistic scenario but allows 
us to plot all power spectra with the same axes. 




0.1 1 10 

Time-lag (ms) 

FIG. 8. Mean-squared-displacement for the Ornstein- 
Uhlenbeck, optical-trap, and Brownian-motion limits de- 
scribed in Sec. VII A, VII B, and VII C. Simulation settings 
and legends are the same as in Fig. 7, except green/blue lines 
show slopes of one/two. 



there are three dimension-full parameters and two 
dimension-full variables. Dividing through with one of 
the parameters there are only two left. From these 
we can form a characteristic time and a characteristic 
length. By expressing times and lengths in terms of 
these, we have removed all parameters from the equa- 
tion. In other words, there is only a single unique equa- 
tion. 

What is a natural choice for the characteristic time 
and length? By considering the physics, we see that 
switching off the driving-force will result in a solution to 
the dynamics that is an exponentially decaying distance 
from x(t) to x — 0. The characteristic time, j/n, for this 
decay is the characteristic time for the dynamics and it 
makes sense to measure time in units of this time. With 
the driving-force on, a statistical equilibrium establishes 
itself. In this equilibrium the distribution of positions 
has a certain width that is temperature-dependent be- 
cause the driving-force is thermal. This width U (x 2 ), is 
a natural unit in which to measure distances and gives 
us a position variable that is dimensionless. 

The basic idea then, is to express the independent 
variable t, and dependent variable x in terms of di- 
mensionless ones and then use the freedom we have 
to choose parameters to make as many pre-factors as 
possible equal to unity [4]. Above we used a physi- 
cal approach. In a more mathematical approach the 
recipe is to write x = x c \ and t = t c 9 where x c and t c 
arc dimension-full "characteristic" quantities, whereas 
X and 9 are the new dimensionless dependent and in- 
dependent variables respectively. Dividing by 7 we first 
get 

— -JK + ~ X oX = "/therm (0) (87) 

t c d0 7 7 

where 

/therm(0) = F thcrm (t c 9) = y/2k B Tj/t c T}(9) (88) 

with 

f,{6) = Vt~ cV (t c e) (89) 

and 

mm')) = He -n m = o. m 

Note, that rj is a dimensionless generalized function of 
the dimensionless parameter 9, whereas 77 is a general- 
ized function of t and has dimension Next, by 
choosing 

t c = 7/K and x c — \j2ksT / k (91) 

we get 

^+X = V(9), (92) 

as our dimensionless dynamic equation. Without hav- 
ing actually solved the original equation we have gained 
three insights: (i) The system has a characteristic time- 
scale of j/k; (ii) the system has a characteristic length- 
scale of y/2k^T/ n (compare to the MSD); and (iii) there 
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is only one OT-equation, not a whole family for different 
values of 7 and k. 

Since the OU-process for the velocity 



mi) + 7U = F t herm(i); V = X 



(93) 



is mathematically identical to the OT process for the 
positions we directly have the dimensionless pendant 



dv 



(94) 



by writing v = v c v, t = t c 9, and choosing 



t c = m/7 and v c = \j2k-QT/m . 



(95) 



For completeness, we also give the dimensionless ver- 
sion of the original dynamics. The dimension- full equa- 
tion 



TUX + JX + KX = F t herm(*) 



(96) 



can be rewritten as the universal oscillator equaition 

(97) 



where A = 7/ ' \fmk is the damping ratio, x c = 
2kBTj/V niK 3 , and t 2 c = m/n. That is, there are 
three dimension-full parameters and a family of solu- 
tions parametrized by A. 

To arrive at dimensionless versions of the discrete 
equations we can proceed to solve the above dimen- 
sionless equations as before. Alternatively, we can sim- 
ply take the already known dimension-full results and 
nondimcnsionalizc them. 



IX. SUMMARY AND CONCLUSIONS 

We examined the dampened harmonic oscillator and 
three of its physical limits: The mass-less case (opti- 
cal trap), the free case (the Einstein-Ornstein-Uhlcnbeck 
theory for Brownian motion) , and the mass-less free case 
(Einstein's Brownian motion). By solving the system's 
dynamical equations for an arbitrary time-lapse At, ex- 
act analytical expressions were derived for the changes 
in position and velocity during such a time-lapse. With 
these expressions, exact simulations of the dynamics are 
then possible — with an accuracy that is independent of 
the duration of the time-lapse. In contrast, a numerical 
simulation, using Euler-integration or similar schemes, 
is exact only to first or second order in At [18]. 

We gave exact analytical expressions for power- 
spectral forms, mean-squared-displacements, and 
correlation-functions that can be fitted (see [19] before 
undertaking a lcast-squares-fit) to data obtained from 
time-lapse recording of a system with dynamics similar 
to the dampened harmonic oscillator or one of its three 
physical limits described here. The effect of finite 
sampling rates (aliasing) were also discussed. 

The effect on the power spectrum of velocity estima- 
tion from position-data (secant-velocity) was treated for 
the case of free diffusion of a massive particle. Approxi- 
mate as well as exact corrective factors and expressions 
were given. Finally, we discussed some of the advan- 
tages of dimensional analysis and of nondimensionaliza- 
tion of the dynamic equations. Throughout, we pointed 
out when power spectral analysis makes sense (bounded 
process) or not, which of the statistical measures that 
depend on the sampling-frequency, and which may be 
described by the simpler continuous-time theory. 
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